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Abstract. We present a numerical study of the hydrodynamics in the final stages of 
inspiral in a black hole— neutron star binary, when the separation becomes comparable 
to the stellar radius. We use a Newtonian three-dimensional Smooth Particle 
Hydrodynamics (SPH) code, and model the neutron star with a soft (r = 5/3) 
polytropic equation of state and the black hole as a Newtonian point mass which 
accretes matter via an absorbing boundary at the Schwarzschild radius. Our initial 
conditions correspond to tidally locked binaries in equilibrium. The dynamical 
evolution is followed for approximately 23 ms, and in every case for T = 5/3 we 
find that the neutron star is tidally disrupted on a dynamical timescale, forming a 
dense torus around the black hole that contains a few tenths of a solar mass. A 
nearly baryon— free axis is present in the system throughout the coalesence, a fireball 
expanding along that axis would avoid excessive baryon contamination and could 
give rise to a modestly beamed gamma-ray burst. 



Introduction 



Theoretical studies of the binary coalescence of a neutron star with a black hole are 
as venerable as the Texas Symposia (Wheeler 1971 0). Such coalescences have been 
implicated in the creation of heavy nuclei through the r-process and were suspected of 
causing the violent outbursts known as gamma-ray bursts (Lattimer & Schramm 1974, 
1976 [|[ [|; they are also candidate sources for the gravitational radiation detectors 
currently coming of age, as the double neutron star binaries have been for a long time 
(Clark & Eardley 1977 §). 

No binary system composed of a black hole and a neutron star has as yet been 
identified as such by the astronomers. However, theoretical estimates and studies of 
binary stellar evolution predict that such systems are created, and that the expected 
rate of their coalescence may be between about one per one hundred thousand and one 
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per a million years per galaxy (Lattimer & Schramm 1974, 1976 Q, [|; Narayan, Piran 
& Shemi 1991 §); Tutukov & Yungelson 1993 §]; Lipunov, Postnov & Prokhorov 

1997 0). This rate is about right to explain the r-process and gamma-ray bursts 
(although it remains to be demonstrated that the phenomena are related to the 
coalescence), and would give a very satisfactory event rate for future gravity wave 
detectors (10 to 100 per year out to 200 Mpc). 

The expected outcome of the coalescence ranged from the neutron star "plunging" 
into the black hole (Bardeen, Press & Teukolsky 1972 ||), through the disruption of 
the neutron star and the formation of a transient accretion stream (Wheeler 1971 [1], 
Lattimer & Schramm 1976 j|) or of an accretion disk/torus (Paczyhski 1991 [9], 
Jaroszyhski 1993 p0| , Witt et al. 1994 ]llj| ), to rapid growth of the binary separation 
in a steady Roche-lobe overflow scenario (Blinnikov et al. 1984 p2[ , Portegies Zwart 

1998 [fL3)). It seemed important to investigate the hydrodynamics of the process 
numerically, to determine which, if any, of these outcomes are likely. 

We have initially performed Newtonian simulations treating the neutron star as 
a stiff polytrope (Lee & Kluzniak 1995 motivated by a desire to determine 

the timescale of the coalescence with the black hole and to investigate the spatial 
distribution of the matter lost by the star. These questions were of relevance to the 
theory of gamma-ray bursts, which seems to allow the coalescence to power the burst 
in the accepted relativistic shock model, provided that the environment into which 
the fireball expands has a very small number of baryons at least in some directions 
(Meszaros & Rees 1992, 199 3 [l5[ and that the explosive event is highly variable 
in time (Sari & Piran 1997 Jl7[). We found the coalescence of a neutron star and a 
black hole promising in these respects, at least in our Newtonian model (Kluzniak & 
Lee 1998 @). 

However, the stiff polytrope (L = 3), used by us in the work cited above, differed 
in one crucial respect from a neutron star — it responded to mass loss by shrinking, 
instead of expanding. For this reason we have repeated the simulation with a softer 
polytrope, and we report the results below, and elsewhere (Lee & Kluzniak 1999 |Ts|). 



2. Numerical method 



For the simulations presented here, we have used the method known as Smooth 
Particle Hydrodynamics. Our code is three-dimensional and is completely Newtonian, 
although removal of angular momentum by gravitational radiation was included, as 
described below. This method has been described often, we refer the reader to 
Monaghan (1992) @ for a review of the principles of SPH, and to Lee (1998) [fl| 
and Lee & Kluzniak (1998) for a detailed description of our own code. 

Here, we model the neutron star via a polytropic equation of state, P = Kp T with 
r = 5/3. The unperturbed (spherical) neutron star has a radius R=13.4 km and mass 
M=1.4M . The black hole (of mass Mbh) is modeled as a Newtonian point mass, 
with a potential $bh(?*) = —GMss/r. We model accretion onto the black hole by 
placing an absorbing boundary at the Schwarzschild radius (rsch — 2GA/bh/c 2 ). Any 
particle that crosses this boundary is absorbed by the black hole and removed from 
the simulation. The mass and position of the black hole are continously adjusted so 
as to conserve total mass and total momentum. 

Initial conditions corresponding to tidally locked binaries in equilibrium are 
constructed in the co-rotating frame of the binary for a range of separations r and 
a given value of the mass ratio q — Mns/Mbh (Rasio & Shapiro 1994 ]23] ]; Lee 
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& Kluzniak 1998 p2|]). The neutron star is modeled with TV = 17,256 particles 
at the start of the calculation in every case presented here. We have calculated 
the gravitational radiation signal emitted during the coalescence in the quadrupole 
approximation, and refer the reader to Lee & Kluzniak (1999) |l9| for details. 

We have included a term in the equations of motion that simulates the effect 
of gravitational radiation reaction on the components of the binary system in the 
quadrupole approximation (see Landau & Lifshitz 1975 This formulation of the 

gravitational radiation reaction has been used in SPH before (Davies et al. 1994 p5| , 
Zhuge et al. 1996 Rosswog et al. 1999 |27|) in the case of merging neutron 

stars, and it is usually switched off once the stars come into contact, when the point- 
mass approximation clearly breaks down. We are assuming then, that the polytrope 
representing the neutron star can be considered as a point mass for the purposes of 
including radiation reaction. Clearly, the validity of this assumption must be verified 
a posteriori when the simulation has run its course. If the neutron star is disrupted 
during the encounter with the black hole, this radiation reaction must be turned off, 
since our formula would no longer give meaningful results. We have adopted a switch 
for this purpose, as follows: if the center of mass of the SPH particles comes within 
a prescribed distance of the black hole (effectively a tidal disruption radius), then the 
radiation reaction is turned off. This distance is set to r t idai = CR^bh/Mns) 1 ^ 3 > 
where C is a constant of order unity. 

3. Results 

3.1. Evolution of the binary 

To allow comparisons of results for differing equations of state, we have run simulations 
with the same initial binary mass ratios as previously explored for T = 3 (Lee & 
Kluzniak 1998 (2^]), namely q—1, q—0.8 and g=0.31. Additionally we have examined 
the case with mass ratio q=0.1. Equilibrium sequences of tidally locked binaries 
were constructed for a range of initial separations, terminating at the point where 
the neutron star overflows its Roche Lobe (at r — trl). In Figure |l|a we show the 
variation of total angular momentum J for one of these sequences as a function of 
binary separation (solid line). Following Lai, Rasio & Shapiro (1993b) we have 
also plotted the variation in J that results from approximating the neutron star as 
compressible tri-axial ellipsoid (dashed lines) and as a rigid sphere (dotted lines). In 
all cases, the SPH results are very close to the ellipsoidal approximation until the 
point of Roche-Lobe overflow. This result is easy to understand if one considers that 
the softer the equation of state, the more centrally condensed the neutron star is and 
the less susceptible to tidal deformations arising from the presence of the black hole. 

For r = 3 (Lee & Kluzniak 1998 |22|]), the variation in angular momentum as a 
function of binary separation was qualitatively different (for high mass ratios) from our 
present findings for T = 5/3. For q=l and q—0.8, total angular momentum attained 
a minimum at some critical separation before Roche-Lobe overflow occurred. This 
minimum indicated the presence of a dynamical instability, which made the binary 
decay on an orbital timescale. This purely Newtonian effect arose from the tidal 
interactions in the system (Lai, Rasio & Shapiro 1993a Q). In the present study, we 
expect all orbits with initial separations r > r^L to be dynamically stable. 

For polytropes, the mass-radius relationship is R cx A/( r ~ 2 )/( 3r ~ 4 ). For T=5/3, 
this becomes R cx A/ -1 / 3 . Thus, the polytrope considered here responds to 
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Figure 1 . (a) Total angular momentum J as a function of binary separation r along 
the equilibrium sequence for a mass ratio g=0.31 (run D). The solid line is the result 
of the SPH calculation, the dashed line results from approximating the neutron star 
as a compressible tri-axial ellipsoid and the dotted line from approximating it as a 
rigid sphere. (b) Binary separation during the dynamical runs for various initial mass 
ratios (see Table ph. The vertical line across each curve (if present), indicates the 
time when gravitational radiation reaction was switched off. 



mass loss by expanding, as do neutron stars modeled with realistic equations of 
state (Arnett & Bowers 1977 j}0))-the dynamical disruption of the star reported 
below seems to be related to this effect. For the polytropic index considered in Lee & 
Kluzniak (1998) [p2j, the star was not disrupted (see also Lee & Kluzniak 1995 @; 
1997 Kluzniak & Lee 1998 @), but we find no evidence in any of our dynamical 
calculations for a steady mass transfer in the binary, such as the one suggested in the 
literature (e.g. Blinnikov et al. 1984 @; Portegies Zwart 1998 @). 

Using the quadrupole approximation, one can compute the binary separation as a 
function of time for a point-mass binary, and obtain 

r = n (1 - t/t f 4 , (1) 

with 1 = 256G 3 MbhM ns (M B h + A/ N s)/(5r l 4 c 5 ). Here n is the separation at t=0. 
For black hole-neutron star binaries studied here, the timescale for orbital decay 
because of angular momentum loss to gravitational radiation, in, is on the order of 
the orbital period, P (for q=l, at an initial separation ri=2.7R we find io=6.5 ms and 
P=2.24 ms). 

3.2. Run parameters 

In Table [l] we present the parameters distinguishing each dynamical run we performed. 
All times are in milliseconds and all distances in kilometers. The runs are labeled 
with decreasing mass ratio (increasing black hole mass), from q=l down to q=0.1. 
All simulations were run for the same length of time, tti na i = 22.9 ms (this covers 
on the order of ten initial orbital periods for the mass ratios considered). The initial 
separation for each dynamical run is given as r,-, and the separation at which Roche 
Lobe overflow from the neutron star onto the black hole occurs is given by r^L- 
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Table 1. Important parameters for each run 



Run 


q 


r RL (km) 


ri(km) 


t rad (ms) 


M ms ) 


N 


A 


1.00 


35.78 


36.18 


3.49 


22.9 


17,256 


B 


0.80 


37.65 


38.19 


3.43 


22.9 


17,256 


D 


0.31 


48.11 


48.24 


1.83 


22.9 


17,256 


E 


0.10 


67.13 


67.67 


1.09 


22.9 


17,256 



The fifth column in Table |l| shows the value of t ra d, when radiation reaction is 
switched off according to the criterion established in section 0. We note here that run 
E is probably at the limit of what should be inferred from a Newtonian treatment 
of such a binary system. The black hole is very large compared to the neutron star, 
and the initial separation (r^ = 67.87 km) is such that the neutron star is well within 
the innermost stable circular orbit around a Schwarzschild black hole of the mass 
considered. 



3.3. Morphology of the mergers 

The initial configurations are close to Roche Lobe overflow, and mass transfer from 
the neutron star onto the black hole starts within one orbital period for all runs 
presented here. In every run the binary separation (solid lines in in Figure [j]) 
initially decreases due to gravitational radiation reaction. For high mass ratios, (runs 
A, B) the separation decays faster than what would be expected of a point-mass 
binary. This is also the case for a stiff equation of state, in black hole-neutron star 
mergers (Lee & Kluzniak 1998 p2[) as well as in binary neutron star mergers (Rasio 
& Shapiro 1994 ^3|), and merely reflects the fact that hydrodynamical effects are 
playing an important role. For the soft equation of state studied here, there is the 
added effect of 'runaway' mass transfer because of the mass-radius relationship (see 
section 3.1). For run C, the solid and dashed lines in Figure [l]b follow each other very 



closely, indicating that the orbital decay is primarily driven by angular momentum 
losses to gravitational radiation. For run E, the orbit decays more slowly than what 
one would expect for a point-mass binary. This is explained by the fact that there is a 
large amount of mass transfer (10% of the initial neutron star mass has been accreted 
by t — t r ad in this case) in the very early stages of the simulation, substantially altering 
the mass ratio in the system (the dashed curves in Figure 0b are computed for fixed 
masses; at constant total mass, note that from equation |], lowering the mass ratio in 
the system slows the orbital decay for q < 0.5). 

The general behavior of the system is qualitatively similar for every run. Figure || 
shows density contours in the orbital plane (left panels) and in the meridional plane 
containing the black hole (right panels) for run D at t = 5.73 ms and t — tj — 22.9 ms. 
The neutron star becomes initially elongated along the binary axis and an accretion 
stream forms, transferring mass to the black hole through the inner Lagrange point. 
The neutron star responds to mass loss and tidal forces by expanding, and is tidally 
disrupted. An accretion torus forms around the black hole as the initial accretion 
stream winds around it. A long tidal tail is formed as the material furthest from the 
black hole is stripped from the star. Most of the mass transfer occurs in the first 
two orbital periods and peak accretion rates reach values between 0.5 M /ms and 
1.2MQ/ms (see Figure ||). The mass accretion rate then drops and the disk becomes 
more and more azimuthally symmetric, reaching a quasi-steady state by the end of 
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Figure 2. Density contours after the disruption of the polytrope for q = 0.31 
(Run D) in the equatorial plane (left panels), and in the meridional plane containing 
the black hole (right panels) at t=5.73 ms (top panels) and 4=22.9 ms (bottom 
panels). The axes are labeled in units of the initial (unperturbed) neutron star 
radius R=13.4 km. All contours are logarithmic and equally spaced every 0.5 dex. 
The lowest contour is at logp/po = — 6 (po = 1.14 X 10 15 g cm -3 ), and bold contours 
are plotted at logp/po = —6, —5, —4, —3. 



the simulations. 

We show in Figure | the various energies of the system (kinetic, internal, 
gravitational potential and total) for run D. The dramatic drop in total internal energy 
reflects the intense mass accretion that takes place within the first couple of orbits. 
Figure ^ also shows [panel (b)] the total angular momentum of the system for runs A, 
B, D and E (the only contribution to the total angular momentum not plotted is the 
spin angular momentum of the black hole, see below) . Angular momentum decreases 
for two reasons. First, if gravitational radiation reaction is still acting on the system, 
it will decrease approximately according to the quadrupole formula. Second, whenever 
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Figure 4. (a) Various energies of the system as a function of time for run D. 
The kinetic (T), internal (U), gravitational potential (W) and total (E) energies are 
indicated, (b) Angular momentum as a function of time for every run. 



matter is accreted by the black hole, the corresponding angular momentum is removed 
from our system. In reality, the angular momentum of the accreted fluid would increase 
the spin of the black hole. We keep track of this accreted angular momentum and 
exhibit its value in Table || as the Kerr parameter of the black hole. This shows up as 
a decrease in the total value of J. 



3.4- Accretion disk structure 

In Table ^| we show several parameters pertaining to the final accretion structure 
around the black hole for every run. The mass that has been accreted by the black 
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Table 2. Accretion disk structure 



Run 


q 


M dlsk (M Q ) 


M acc (M ) 


M ma ,(M /nis) 


M/w(M /s) 


T dlsk (ms) 


A 


1.00 


0.263 


1.092 


0.831 


4.88 


54.1 


B 


0.80 


0.277 


1.078 


0.733 


6.11 


46.9 


D 


0.31 


0.316 


1.036 


0.549 


4.88 


63.1 


E 


0.10 


0.101 


1.204 


1.160 


2.44 


46.9 



Table 3. 



Run 


Jbhc/GM| h 


0-3 


6-A 


0- 5 


A 


0.517 


20 


12 


8 


B 


0.497 


25 


10 


3 


D 


0.173 


40 


21 


10 


E 


0.114 


52 


42 


32 



Q- n is the half-angle (in degrees) of a cone above the black hole and along the rotation 
axis of the binary that contains a mass M = 1.4 x 10~™M Q . 




Figure 5. (a) Azimuthally averaged profiles for the density, p, and the specific 
internal energy, u (n/10 8 is plotted), of the accretion torus at t = tf for run D. The 
inner edge of the curves is at r = 2rg c f l . At this stage in the simulation the torus is 
close to being azimuthally symmetric, (b) Distribution of specific angular momentum 
j for runs A, B, D and E at t = tf (solid lines, A, B, D, E) and the specific angular 
momentum of a Kcplerian accretion disk for the same black hole mass (dashed lines, 
A, B, D, E). 
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hole is denoted by M acc . The disk settles down to a fairly azimuthally symmetric 
structure within a few initial orbital periods (except for the long tidal tail, which 
always persists as a well-defined structure) , and there is a baryon-free axis above and 
below the black hole in every case (see below). We have calculated the mass of the 
remnant disk, Mdisk, by searching for the amount of matter that has sufficient specific 
angular momentum j at the end of the simulation to remain in orbit around the black 
hole (as in Ruffert & Janka 1998 p2[). This material has j > jcrit — \/&GM t /c, 
where M t is the total mass of the system. By the end of the simulations, between 
70% and 80% of the neutron star has been accreted by the black hole. It is interesting 
to note that the final accretion rate (at t = tf) appears to be rather insensitive to 
the initial mass ratio, and is between 2.4 M© s _1 and 6.1 M s _1 . From this final 
accretion rate we have estimated a typical timescale for the evolution of the accretion 
disk, Tdisk — Mdisk/ M final- Despite the difference in the initial mass ratios and the 
typical sizes of the disks, the similar disk masses and final accretion rates make the 
lifetimes comparable for every run. 

We have plotted azimuthally averaged density and internal energy profiles in 
Figure ^ for run D. The specific internal energy is greater towards the center of 
the disk, and flattens out at a distance from the black hole roughly corresponding 
the density maximum, at u ~ 3 x 10 18 erg g _1 , or 3.1 MeV/nucleon, and is largely 
independent of the initial mass ratio. The inner regions of the disks have specific 
internal energies that are greater by approximately one order of magnitude. 

Additionally panel (b) in the same figure shows the azimuthally averaged 
distribution of specific angular momentum j in the orbital plane for all runs. The curves 
terminate at ri n = Irgch- Pressure support in the inner regions of the accretion disks 
makes the rotation curves sub-Keplerian, while the flattening of distribution marks 
the outer edge of the disk and the presence of the long tidal tail (see Figure ||) , which 
has practically constant specific angular momentum. 

The Kerr parameter of the black hole, given by a = Jbhc/GM| h , is also shown 
in Table [|. We have calculated it from the amount of angular momentum lost by the 
fluid via accretion onto the black hole (see Figure ||b), assuming that the black hole 
is not rotating at t — 0. The specific angular momentum of the black hole is smaller 
for lower mass ratios simply because the black hole is initially more massive when q 
is smaller. 

It is of crucial importance for the production of GRBs from such a coalescence 
event that there be a baryon-free axis in the system along which a fireball may 
expand with ultrarelativistic velocities (Meszaros & Rees 1992, 1993 JlB], 0). We 
have calculated the baryon contamination for every run as a function of the half- 
angle A9 of a cone directly above the black hole and along the rotation axis of the 
binary that contains a given amount of mass AM. Table || shows these angles (in 
degrees) for AM/M = 1.4 x 10~ 3 , 1.4 x 10~ 4 , 1.4 x 10~ 5 . There is a greater amount 
of pollution for high mass ratios (the disk is geometrically thicker compared to the 
size of the black hole) , but in all cases only modest angles of collimation are required 
to avoid contamination. We note here that the values for 9-$ are rough estimates at 
this stage since they are at the limit of our numerical resolution in the region directly 
above the black hole. 
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4. Discussion 



The numerical simulations reported here were quasi-newtonian. An important caveat 
to keep in mind is that inclusion of general relativistic effects may lead to results 
qualitatively different from even post-newtonian treatment (Wilson, Mathews & 
Marronetti 1996 @). 

Our results indicate that the outcome of the binary coalescence depends on the 
nature of the star orbiting the black hole. As reported previously (Lee & Kluzniak 
1995 @, 1998 Kluzniak & Lee 1998 ]!§]), when we modeled the star as a 

polytrope with adiabatic index T = 3, the coalescence appeared to be an intermittent 
process in which the core of the polytrope survives the initial encounters and increases 
its separation from the black hole, thus extending the merger to possibly ~ 0.1 s. For 
the softer polytrope discussed here (r = 5/3), the star is disrupted completely in a 
few milliseconds and all that remains after the initial mass transfer is an accretion 
disk, containing no more than 1/5 of the initial mass, and some ejecta. Perhaps the 
current simulation with F = 5/3, is the more realistic one, because for this polytrope 
dM/dR < 0, as for physical models of neutron stars (e.g. Arnett & Bowers 1977 

In agreement with earlier suggestions (Lattimer & Schramm 1974, 1976 fi Ej|), we 
have found that some matter will be ejected from the system, in an amount sufficient 
to account for the abundance of the r-process nuclei (assuming the r-process does 
indeed occur during the merger). 

The binary coalescence of a neutron star with a black hole remains an attractive 
theoretical source of gamma-ray bursts. The energy requirements for at least one 
recently observed burst are so severe, if emission is isotropic (e.g. Kulkarni et al. 
1998 1 34 ) , that some degree of beaming seems desirable. According to the simulations 
presented here, ultrarelativistic flows are possible in the post-merger system only along 
the rotational axis of the system in a solid angle of about 0.1 steradian. Proper 
inclusion of neutrino transport may change this angle somewhat. The rather short 
(~ 50 ms) accretion timescale of the remnant disk reported, does not include possible 
interaction between the disk and the black hole (Blandford & Znajek 1977 p5|). 
In fact, the appearance in the simulation of a substantial toroidal disk around the 
black hole is encouraging, as it may allow the black hole spin to be extracted by the 
Blandford-Znajek mechanism, possibly powering in this manner the gamma-ray burst 
fireball (Meszaros & Rees 1997 
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